home *** CD-ROM | disk | FTP | other *** search
/ Aminet 1 (Walnut Creek) / Aminet - June 1993 [Walnut Creek].iso / aminet / misc / sci / ephem_src_4_28.lha / sun.c < prev    next >
C/C++ Source or Header  |  1992-04-17  |  2KB  |  56 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include "astro.h"
  4.  
  5. /* given the modified JD, mjd, return the true geocentric ecliptic longitude
  6.  *   of the sun for the mean equinox of the date, *lsn, in radians, and the
  7.  *   sun-earth distance, *rsn, in AU. (the true ecliptic latitude is never more
  8.  *   than 1.2 arc seconds and so may be taken to be a constant 0.)
  9.  * if the APPARENT ecliptic longitude is required, correct the longitude for
  10.  *   nutation to the true equinox of date and for aberration (light travel time,
  11.  *   approximately  -9.27e7/186000/(3600*24*365)*2*pi = -9.93e-5 radians).
  12.  */
  13. sunpos (mjd, lsn, rsn)
  14. double mjd;
  15. double *lsn, *rsn;
  16. {
  17.     double t, t2;
  18.     double ls, ms;    /* mean longitude and mean anomoay */
  19.     double s, nu, ea; /* eccentricity, true anomaly, eccentric anomaly */
  20.     double a, b, a1, b1, c1, d1, e1, h1, dl, dr;
  21.  
  22.     t = mjd/36525.;
  23.     t2 = t*t;
  24.     a = 100.0021359*t;
  25.     b = 360.*(a-(long)a);
  26.     ls = 279.69668+.0003025*t2+b;
  27.     a = 99.99736042000039*t;
  28.     b = 360*(a-(long)a);
  29.     ms = 358.47583-(.00015+.0000033*t)*t2+b;
  30.     s = .016751-.0000418*t-1.26e-07*t2;
  31.     anomaly (degrad(ms), s, &nu, &ea);
  32.     a = 62.55209472000015*t;
  33.     b = 360*(a-(long)a);
  34.     a1 = degrad(153.23+b);
  35.     a = 125.1041894*t;
  36.     b = 360*(a-(long)a);
  37.     b1 = degrad(216.57+b);
  38.     a = 91.56766028*t;
  39.     b = 360*(a-(long)a);
  40.     c1 = degrad(312.69+b);
  41.     a = 1236.853095*t;
  42.     b = 360*(a-(long)a);
  43.     d1 = degrad(350.74-.00144*t2+b);
  44.     e1 = degrad(231.19+20.2*t);
  45.     a = 183.1353208*t;
  46.     b = 360*(a-(long)a);
  47.     h1 = degrad(353.4+b);
  48.     dl = .00134*cos(a1)+.00154*cos(b1)+.002*cos(c1)+.00179*sin(d1)+
  49.                                 .00178*sin(e1);
  50.     dr = 5.43e-06*sin(a1)+1.575e-05*sin(b1)+1.627e-05*sin(c1)+
  51.                         3.076e-05*cos(d1)+9.27e-06*sin(h1);
  52.     *lsn = nu+degrad(ls-ms+dl);
  53.     *rsn = 1.0000002*(1-s*cos(ea))+dr;
  54.     range (lsn, 2*PI);
  55. }
  56.